library(lme4)
#> Loading required package: Matrix
library(glmmTMB)
library(dplyr) ## for mutate_at, %>%
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(performance)

#Build example data
x <- c("A", "B", "C", "D")
time <- rep(x, each=20, times=3) #time factor
y <- c("exposed", "ref1", "ref2")
lake <- rep (y, each=80)  #lake factor
set.seed(123)
min <- runif(n=240, min=4.5, max=5.5) #mins used in model offset
set.seed(123)
count <- rnbinom(n=240,mu=10,size=100) #randomly generated negative binomial data

#make data frame
dat <- as.data.frame(cbind(time, lake, min, count)) 
dat <- dat |> 
   mutate_at(c('min', 'count'), as.numeric)

#remove one combination of factors to make example rank deficient (all observations from time A and lake ref1)
dat2 <- filter(dat, time!="A" | lake !="ref1")

model <-glmmTMB(count~time*lake, family=nbinom1,
                      control = glmmTMBControl(rank_check = "adjust"),
                      offset=log(min), data=dat2)
#> dropping columns from rank-deficient conditional model: timeD:lakeref1
check_overdispersion(model) |> plot()

set.seed(3)
mu <- rpois(500, lambda = 3)
x <- rnorm(500, mu, mu * 3) |> ceiling() |> pmax(0)

quine.nb1 <- MASS::glm.nb(x ~ mu)
check_overdispersion(quine.nb1) |> plot()
#> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

set.seed(101)
d <- data.frame(x = runif(1000), f = factor(sample(1:200, size = 1000, replace = TRUE))) # modified for more random effects
suppressMessages(
  d$y <- simulate(
    ~ x + (1 | f),
    family = poisson,
    newdata = d,
    newparams = list(theta = 1, beta = c(0, 2))
  )[[1]]
)
m1 <- glmer(y ~ x + (1 | f), data = d, family = poisson)
check_overdispersion(m1) |> plot()

docvisit <- datawizard::data_read("~/../Downloads/docvisit.txt")

mp <- glmmTMB(
  doctorco ~ sex + illness + income + hscore, 
  data = docvisit,
  family = poisson()
)
check_overdispersion(mp) |> plot()

mnb <- glmmTMB(
  doctorco ~ sex + illness + income + hscore,
  data = docvisit,
  family = nbinom2()
)
check_overdispersion(mnb) |> plot()

mzip <- glmmTMB(
  doctorco ~ sex + illness + income + hscore,
  ziformula = ~age,
  data = docvisit,
  family = poisson()
)
check_overdispersion(mzip) |> plot()

mzinb <- glmmTMB(
  doctorco ~ sex + illness + income + hscore,
  ziformula = ~age,
  data = docvisit,
  family = nbinom2()
)
check_overdispersion(mzinb) |> plot()

mzinbd <- glmmTMB(
  doctorco ~ sex + illness + income + hscore + age,
  ziformula = ~ sex + illness + income + hscore + age,
  dispformula = ~ sex + illness + income + hscore + age,
  data = docvisit,
  family = nbinom2()
)
check_overdispersion(mzinbd) |> plot()